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^j' ' Abstract 

t^ , From the numerical point of view, given a set X C K" of s points whose 

coordinates are known with only limited precision, each set X of s points 
whose elements differ from those of X of a quantity less than the data 

^^ • uncertainty can be considered equivalent to X. We present an algorithm 

'^H ' that, given X and a tolerance £ on the data error, computes a set Q of 

ri , polynomials such that j;ach element of Q "almost vanishing" at X and at 

all its equivalent sets X. Even if Q is not, in the general case, a basis of 
the vanishing ideal 2r(X), we show that, differently from the basis of 2r(X) 
that can be greatly influenced by the data uncertainty, Q can determine a 
geometrical configuration simultaneously characterizing the set X and all 

»vj , its equivalent sets X. 

K^ ' Keywords: Vanishing ideal, border and Grobner bases, limited precision data. 

^ : 1 Introduction 

r*~" ' Let P — R[xi, . . . , Xn] be the polynomial ring in n indcterminates over the reals 



u 



C^ 



o 

oo 

o 



c^ 



and let X = {pi, . . . ,ps} be a finite set of points of M". 

It is well known [21 |TU] that the vanishing ideal I(X) C P of all polynomials 
which vanish at X can be described by a Grobner basis [3, if a term ordering 
is chosen, or by a border basis, if an appropriate basis of the quotient space 
^ : P/I(X) is given. 

^ ' However, it is also well known that small perturbations of the points of X 

can cause structural changes in the bases of I(X) [TUl [T3] as illustrated in the 
following example. 

Example 1.1 Let a be the DcgLex term ordering with x > y. Given the set 
of points X = {(1, 1), (3, 2), (5.1, 3)}, the a-Grobner basis GB of the vanishing 
ideal I(X) is given by: 

j/2 - 20a; + 37y - 18 
GB^ ■{ xy^ 43x + Sly - 39 

c2-90.1x+172.2y-83.1 
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The set GB is also the border basis of X(X), founded on the set O = {1, y, x} 
whose residue classes span P/1{X). 

A slightly perturbation of the point (5.1,3) leads to a new set of points 
X = {(1, 1), (3, 2), (5, 3)}. The cr-Grobncr basis GB of the vanishing ideal I{X) 
is completely different from GB: 



^^ — "> ,3 «„,2 



x-2y+l, 



Further, all the border bases of 2r(X) also present a structural discontinuity, 
since the residue classes of the set O do not span the space P/I{X). <0> 

In the previous example the structural changes happen since the input points X 
are almost aligned, while the slightly perturbed points X are exactly aligned. 
This example also illustrates that, if we deal with a set X of points known with 
limited precision, the exact bases of X(X) could not highlight some pleasant 
geometrical properties almost satisfied by the points X. 

In this paper we present an algorithm that computes, given a set of points 
known with limited precision, a set of polynomials allowing to recognize if such 
points almost lye on a particularly simple geometrical configuration. 

Given a set X of points whose coordinates are known with limited precision, 
each p of X represents a "cloud" of points: every point p which differs from p 
by less than the data uncertainty can be considered computationally equivalent 
to p. Analogously, an input set obtained from X replacing some p by its per- 
turbation p can be considered an admissible perturbation computationally 
equivalent to X. It is then clear that the knowledge of X with limited preci- 
sion, combined with the structural discontinuity of a basis, points out that a 
significant characterization of I(X) can be a very tricky problem. In fact the 
structure of a basis can drastically change choosing different admissible input 
sets and moreover a blindly choice of a basis can hidden significant geometrical 
properties of X. For this reason exact methods applied to limited precision data 
can produce meaningless results. 

The problem of the characterization of the vanishing ideal of a set of per- 
turbed points has been studied by several authors from different points of view. 
In [llj . Sauer describes a method, suitable for numerical computations, which 
computes a small degree algebraic variety containing the input points. In [7], 
Heldt et al. present an algorithm, based on the singular value decomposition 
of matrices, that computes, without using explicitly the estimation of the data 
error, a set of polynomials which assume particularly small values at the input 
points. 

In [2], Abbott et al. present an algorithm that computes, explicitly using 
the tolerance on the data error, a monomial set O which, in most cases, is a 
basis of P/I(X) and of P/X(X) for all the admissible perturbations X so that 
the O-border basis of all the vanishing ideals X(X) can be obtained. 



Given a set X of limited precision points and a tolerance e on the data error, 
we focus our attention on the possibility of simultaneously characterizing, with 
a set of polynomials, the set X together with all its admissible perturbations. 
To this aim, we present an algorithm that computes an order ideal O and 
a polynomial set tj, whose supports are defined by O, having the following 
properties. 

1. The elements of Q are almost vanishing, w.r.t. the norm of their coefficient 
vectors, at X and at each admissible perturbation X. 

2. For each admissible perturbation X, the set {r(X)|r e O} consists of 
independent vectors, up to the first order error analysis. 

3. For each leading term t oi g ^Q there could be an admissible perturbation 
Xg such that t(Xg) depends on {r(Xg)|r e O}. 

Condition 3 implies that for each g (z G there could exist a polynomial g with 
the same support of g and similar coefficients which vanishes at X^. If it is 
the case, the algorithm determines a geometrical structure, given by g, almost 
satisfied by all the admissible perturbations of X and similar to a geometrical 
structure, given by g, exactly satisfied by the admissible perturbation Xg. 

As illustrated in the numerical examples in Section 6, it can happen that 
there exists a single admissible perturbation X satisfying the previous property 
for all the polynomials g (z G- In this case it is very natural to consider X as a 
possible exact input set, that is the input in absence of data error. Moreover, 
even if in general G is not a basis of X(X), it is analogously natural to consider 
5 as a common characterization of all the admissible perturbations of X. 

Finally, once again as we will show in Section 6, it can happen that X turns 
out to be the exact zero set of the polynomials of G so that, in this case, t/ is a 
Grobner basis of X(X). 

There is an evident open problem regarding the algorithm and its results: 
the existence and possibly the determination of the admissible perturbation X. 
We will show, once again in Section 6, that there are cases when X does not exist. 
Then an open problem is to find conditions for the existence of X and, in case 
of existence, to determine it explicitly. The numerical tests suggest that in case 
of non existence, this can be due to two possible causes: the algorithm does not 
recognize a possible element of O or it detects some geometrical configurations, 
close to the points of X, which are incompatible with each other. The study of 
such open problem will be the subject of our future work. 

The paper improves and formalizes a 2005 Preprint of the author [5] and it 
is organized as follows. Section 2 shows some basic concepts. Section 3 contains 
the description of the algorithm. Section 4 and Section 5 describe the numerical 
properties of O and G, respectively. Finally, Section 6 presents some examples 
illustrating the behaviour of our method. 



2 Preliminaries 

In order to formalize the idea of perturbed points, we recall the definitions of 
empirical point and of admissible perturbation |14[ [2] . 



Definition 2.1 Let p = (ci, . . . ,c„) be a point o/R" and let e — (ei, . . . ,£«), 
with each Si G M+, be the vector of the componentwise tolerances. An empirical 
point p^ is the pair {p,e), where we call p the specified value and e the tol- 
erance. A points p — (ci, . . . ,c„) G M" is called an admissible perturbation 
ofp ifci = Ci + ei, \ei\ < Ei, i = 1, . . . ,n. 

Given a finite set X^ of empirical points all sharing the same tolerance £, we 
can formalize the concept of a set X "equivalent" to X w.r.t. the data accuracy. 

Definition 2.2 Let X*^ = {pi, . . . ,pl} be a set of empirical points with uniform 
tolerance s and with X C K". A set of points X = {pi, . . . ,p^} C R" is called an 
admissible perturbation of X if each pi is an admissible perturbation of pi . 

Finally, we recall (see [Sllin]) some basic concepts related to the polynomial 
ring P = ]R[a:i, . . . ,Xn]. 

Definition 2.3 Let X = {pi, . . . ,ps} be a non-empty finite set of points o/ M" 
and let G — {gi, . . . , gk} be a non-empty finite set of polynomials. 

• The W-linear map evalx : P ^M.^ defined by evalx(/) — (/(pi), • ■ • , fiPs)) 
is called the evaluation map associated to X. For brevity, we write /(X) 
to mean evalx (/)• 

• The evaluation matrix (or vector if k = \) of G associated to X, 
written Mq(K) (or gi(K)), is defined as having entry {i,j) equal to gj{pi). 

Definition 2.4 Let T" be the monoid of power products of P and let O be a 
non-empty subset off". 

• The set O is called an order ideal if O = O, where O is the set of all 
power products in T" which divide some power product of O . 

• Given an order ideal O, the corner set of O is the set 

C[0] ={teT" : t (^ 0,x,\t ^ t/x^ e 0,i ^ 1 . . .n} 

Later on we suppose the reader familiar with the concepts of Grobner basis 
and border basis of a vanishing ideal. Regarding these arguments, the reader is 
referred to the literature (see, among the others, [51 [^[TU]). 



3 The Numerical Algorithm 

Before processing a set X of limited precision points, it is possible to mitigate 
some negative effects of the data uncertainty, replacing with a single represen- 
tative point the elements of X which differ each other by less than the data 
accuracy, since they can be regarded as different perturbations of the same em- 
pirical value. Later on we suppose w.l.o.g. that the set X does not present such 
"redundancy" . If it is not the case, it is possible to preprocess the input data to 
obtain well-separated points, using for instance the algorithms described in [1] 
and included in CoCoALib ^ . Nevertheless the preprocessing of the input data 
is not sufficient to eliminate the instabilities of the exact bases of the vanishing 
ideal I(X), as illustrated in Example ll.il where the points X are well separated. 

We base the construction of our algorithm on the Buchberger-MoUer one [3] 
which computes, given a set X of points and a term ordering cr, the cr-Grobner 
basis GB of X(X) as follows. At each step, if O = {ii, . . . , tk} is the order ideal 
computed at the previous steps, a power product t >„ U is chosen. If the vector 
i(X) is linearly independent of the vectors {ti(X), . . . , ifc(X)}, t is added to O. 

u 

Otherwise, the polynomial g = t— z2i=i ^i^i i^ P^^ vaio GB. Nevertheless, since 
the test of linear dependence is crucially affected by even very small variations 
of the input data, when we deal with points known with limited precision, small 
perturbations of the input data may lead to different choices in the Buchberger- 
MoUer algorithm. 

In order to solve this drawback, we present an algorithm which checks the 
linear dependence in a robust way w.r.t. the data uncertainty. Since every 
admissible perturbation X is computationally equivalent to X, the vector t{%) 
can be considered numerically dependent on {ii(X), . . . , ffe(X)} if there exists an 
X such that i(X) exactly depends on the vectors {ti(X), . . . , tfc(X)}. Formally 
we have the following definition. 

Definition 3.1 Given a set O = {ti, . . . ,tk} and a power product t, the vec- 
tor i(X) numerically depends on {ii(X), . . . , tfc(X)} if there exists an admis- 
sible perturbation 'Kof'K such that the residual /9(X) o/ the least squares problem 
MoQQa = i(X) is a null vector. 

3.1 Sensitivity of the least squares problem 

In order to detect the numerical linear dependency of a set of evaluation vectors, 
we need some results concerning the sensitivity of the least squares problem 

Mo(X)a = t(X) (1) 

First of all we recall some results, based on the componentwise perturbation 
analysis (see [8]), about the sensitivity of a generic least squares problem. 

Given an /i x fc matrix A, we denote by A^ its pseudoinverse, that is 
A+ = {A* A)^^ A* , and by \A\ the matrix consisting of the absolute values of 
the elements of A; given an ft, x fc matrix B, we assume that \A\ < \B\ means 



that the relation holds componentwise. Moreover, given a real value r] -^ 1 we 
denote by 0{vi™)^ m G N, an ft. x fc matrix W{t]) — {wij{r])) (or a real function 
it h = k — 1) such that, for each (i, j), \wij{r])\/i]"^ is bounded near 0. 

Theorem 3.2 Let A and A + AA both be p x q, p > q. full rank matrices and 
let b and b + Ab be two vectors of MP such that \AA\ < rjE and \Ab\ < rjf , 
T] e M+, 77 ^ 1. Consider the least squares problems 

Ax — b with p = b — Ax and 
{A + AA){;x + Ax) ^ b + Ab with p + Ap = b + Ab- {A + AA){;x + Ax) 

We have that 

Ax = A+{Ab-AAx) + [A^A)-^{AAfp + 0{Tf) 
Ap = {I-AA+){Ab-AAx)~{A+Y{AAYp + 0{rf) 

where I is the p x p identity matrix. 

Proof. It is possible (see [S]) to express the least squares problem in the form 
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and so 



{A + AA) 




p + Ap 

X + Ax 



{A + AA)* 
Taking the difference of the previous equations, we have 



b + Ab 
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Ab- AA{x + Ax) 
-iAAYip + Ap) 



Since 



I-AA^ 

A+ 



(A+Y 
-(A'A)-^ 



is the inverse matrix of 



/ 



A 




we obtain 

Ap 
Ax 



(/ - AA+){Ab - AA{x + Ax)) - (A+)*(AA)*(p + Ap) 
A+{Ab - AA{x + Ax)) + {A*A)-\AAY{p + Ap) 



(2) 



Supposing \AA\ < r/E and \Ab\ < rjf, the absolute values of Ax and Ap satisfy 

|Ap| < i]{\l-AA+\{f + E\x + Ax\) + \A+\'E'\p + Ap\) 
\Ax\ < r]{\A+\if + E\x + Ax\) + \{A*A)-^\E'\p + Ap\) 

so we have that 

AAAx = 0{r]^) and (AA)*Ap = 0(77^) 



and the conclusion follows. <C> 

Since we are interested in the behaviour of the least squares problem ([T]) , we 
present an estimation of the sensitivity of the matrix Mq (X) and of the vector 
t(X) to sHght perturbations of the set X. 

Given the power product t ~ x^^ . . . a:^" and the monomial set O, we denote 
by £m = max{ei, i = 1 . . . n}, by deg(a;fe,t) = /3fc the degree of x^ into t, by 
a^i = deg(xfe, t)xf^ . . . xl"-^ ...xt and by dkO = {dkt : i G O}. 

The following result concerns the sensitivity of the evaluation vector i(X). 

Lemma 3.3 Let t he a power product o/T". Given a set X*^ = {pf , . . . ,Pg} o/ 
empirical points and an admissible perturbation X — {pi, . . . ,ps} o/X, we have 
that the vector At = i(X) — i(X) satisfies 

n 

At = ^^fcafci(x) + o(£2^) 
fc=i 

where Ek ~ Diag{ei^k, ■ ■ ■ ,Gs,k) is a diagonal matrix and Ci^k, \&i,k\ < £k, is a 
perturbation the k-th coordinate of pi. 

Proof. First of all we consider a point p = (ci, . . . , c„) £ K" and an admissible 
perturbation p = (ci, . . . , c„) of p w.r.t the tolerance e. Given t — x]^^ . . . xf^" , 
we have that 

t{p) - t{p) = (ci + ei)^i . . . (c„ + e„)'^" - cf^ . . c^ = 

n n n 

Y.^Aci"-' n cft+0{el,)^Y.^kdkt{p) + 0{sl,) 

k=l h=l,h=jtk /c=l 

then we obtain 

n 

t{p^)'^t{p,) = Ye,^k^kt{p^) + 0{elI) 
fc=i 

and so, since t{pi) — t{pi) is the i-th coordinate of t(X) — i(X), we conclude that 

n 

t{X)-t{X) = J2^kdkt{X)+0{el,) 

fc=i 

The following result concerns the sensitivity of the evaluation matrix Mq (X) . 

Lemma 3.4 Let O be an order ideal. Given a set X*^ = {pi . . . ,ps} of empirical 
points and an admissible perturbation X = {pi, . . . ,ps} o/X, we have that the 
matrix AM ~ Mq (X) — Mq (X) satisfies 



AM = ^^fcMa,o(X) + 0(4,) 



fc=i 



where Ek = Diag{ei^k, • ■ ■ , £s,k) is the diagonal matrix of Lemma\K 



Proof. Since the j-th column of Mo{X) is given by ij(X), tj S O, Lemma 
implies that the j-th column of Ma (X) — AIq (X) is 

n 

i,(X) - i,(X) = J2 Ekdkt,(X) + 0{el,) 
fc=i 

The conclusion follows since dktj{X) is the j-th column of the evaluation matrix 
of the set dkO. <} 



The next theorem, based on Theorem 13.21 Lemma 13.31 and Lemma 13. 4|, 
presents a componentwise estimation of the sensitivity of the problem ([T]) to 
the data perturbations. Further, it shows a componentwise upper bound of 
the absolute value of the residual, when there exists an admissible perturbation 
X such that the perturbed least squares problem MoO^)a = i(X) has a zero 
residual. 

Theorem 3.5 Let X"^ be a set of s empirical points and let X be an admissible 
perturbation o/ X. Let O be an order ideal such that Ma{lC) and Me)(X) are 
full rank matrices. Given the least squares problems 

Mo(X)a = t(X) with residual p{X) = t{X) - Mo(^)a 

and _ _ _ _ _ 

Mc)(X)a = i(X) with residual /3(X) = t(X) — Afc'(X)a 

then the vectors Aa — a — a and Ap = /o(X) — p(X) satisfy 

Ap = {I~Mo{X)M+{X))j:l^^Ek{dkt{X)-MaM^)a) 
- (Af+(X))* (ELi Mi,^oi^)Eu) p{X) + 0{el,) 



(3) 



Aa = M+{X)Y:i=iEk{dkt{X)^Ms,o{^)a) 

+ (M*,(X)Mo(X))-i (ELi Mloi^)Ek) p{X) + 0{el,) 



(4) 



Moreover, if there exists an admissible perturbation 'KoflL such that the residual 
p(X) of the least squares problem Mc)(X)S = f(X) is a zero vector, then the 
residual p(X) satisfies 

n 

|p(X)| < |/ - Mo(X)Af+(X)| ^£fc \dkt{X) - Mo,o(X)a| + 0{el,) 

fe=i 

Proof. Since |AM| < emE, \At\ < enf and Mo(X) and Mo{i) have full 
rank, from Theorem 13.21 we obtain 



Ap = (/-Afo(X)Af(^(X))(At-AMa)-(M+(X))*(AM)*p(X) + 0(e|^) 
Aa = M+(X)(At-AMa) + (M^(X)Afo(X))-i(AAf)*p(X) + 0(e^f) 



and so, from Lemma [?751 and Lemma [3^ Equations ^ and ^ follow. 
Moreover, if p(X) is a zero vector from formula ^ we have 



p(X) = -Ap = (Afc)(X)Af+(X) -I)J2Ek idkt{X) - A%o(X)a) + ©(e^,) 

fc=i 

and if we consider the componentwise absolute value of p(X) we obtain 

n 

|p(X)| < \l-Mo{li)M+{X)\Y,\Ek\\dkt{X)-MaM^)a\+Oieli) 

fc=i 

n 

< \I- Afo(X)M+(X)| ^ sk \dkt{X) - A/fa,o(X)a| + 0(4,) ^ 
fe=i 

3.2 The NBM Algorithm 

Theorem 13.51 shows a sufficient condition for the numerical independency of 
i(X) of the columns of Afc)(X). In fact if the residual p(X) of the least squares 
problem Mo{X)a = t(X) satisfies 

n 

|p(X)| > \I - Mo{X)M+iX)\Y,£k\dkt{X) + M9,o{^)a\ + 0{elj) (5) 

fe=i 

then there are no admissible perturbations X of X such that the residual of the 
least squares problem MoOQa ~ t(X) is a null vector. So from Definition 13.11 
it follows that <(X) is numerically independent of {f (X) : r G O}. In particular, 
this implies that if Afe)(X) is a full rank matrix then [Afe)(X)i(X)] is a full 
rank matrix too, for each admissible perturbation X. By exploiting this idea, 
we develop the Numerical Buchberger Moller algorithm, whose main check is 
based on condition ([5]). In particular, since we assume the tolerance on the 
data error is relatively small, we neglect the errors of order 0(e|,) focusing our 
attention on a first order error analysis of the problem. 



The Numerical Buchberger Moller (NBM) Algorithm. 
Input. A set X^ of s empirical points and a term ordering a. 
Output. An order ideal O and a polynomial set Q. 

At the first step O = {1} and G is an empty set. A generic step can be 
described as follows. Let O = {ti, . . . , tk} be the order ideal computed at the 
previous steps and let t be the current power product, t >„ ^i, • ■ • ,ifc, chosen 
according to the strategy of the Buchberger-MoUer algorithm. 

1. Solve the least square problem A/c)(X)a = t(X) and compute the residual 

p(X) =t(X)-A/o(X)a. 

2. If p(X) satisfies 






|p(X)| > |/ - Afo(X)Af+(X)| > efc \dkt{X) - A%o(X)a| (6) 



fe=i 



then put the term t mto the set O. 
3. Otherwise, put the polynomial g = t — X]i=i ^i^i i^^^o Q. ^ 

The NBM algorithm stops after finitely many steps and computes an order 
ideal O, since the strategy to choose the power products to analyze is the same 
as in the Buchberger-MoUer algorithm. 

Note that the term ordering a is only a computational tool for obtaining a 
set O closed under taking divisors. In fact in the general case O is different from 
Oa, the quotient basis determined by the cr-Grobner basis of X(X). Moreover 
it can happen that, for each possible term ordering r, O does not coincide to 
any Or corresponding to the r-Grobner basis of X(X) (see Example 16.41) . For 
this reason any different strategy for building an order ideal can be used in the 
NBM algorithm instead to fix a term ordering. 

4 Properties of the order ideal O 

A first important property of the order ideal O computed by the NBM algorithm 
is its invariance w.r.t. the scaling and the translation of the points X, as shown 
in the following theorem. 

Theorem 4.1 Let X"^ be a set of empirical points with 

X = {pi,...,Ps} Pj = (ci,i,...,Ci^„) and e = (ei, . . . , e„) 

Let X| be the set of scaled empirical points such that X5 = {pi, . . . ,p^} and 

Pj = (diQ^i, . . . ,d„Ci,„) with (di, . . . ,d„) e R" and 5 = (|di|ei, . . . , |d„|e„) 

Let X^ be the set of translated empirical points such that ILt =^ {pi, • ■ • ,Ps} and 

Pi ^ (ci^i +vi,. . . ,Ci^n + Vn) witk (wi ,...,«„)£ K" and T = £ 

Then the NBM algorithm computes the same order ideal O for all the input sets 
X^ X| andX^. 

Proof. We prove that the NBM algorithm computes the same order ideals at 
each step independently of the input sets X*^, X| or X^-. 

At the first step it is true, since O = {1}. Let us suppose that, at the current 
step, with all the three input sets the same order ideal O — {ti, . . . ,tk} has been 
computed and that the term t has to be processed. 

Let us consider the set X5 of the scaled points. 

Given a term r — x^^ . . . x(^", denoting by r{d) ^ dy' . . . d^", we have 

r[d) 
"riPi) = r{d)r{pi) and dkr{pj) = —j—dkr{pi) 

Uk 



10 



so that, denoting by Dq the diagonal matrix Diag{r{d) : r E O), 

i(Xs) = t{d)t{X) and MoO^s) = Mo{X)Do 

dkt{Xs) = ^9fet(X) and MqM^s) = -J-Afa.o(X)i?c) 
rffc rffe 

The least squares problems AIc>{X)a — t{X) and Mc)(X5)a5 = t(Xs) solved 
with input sets X^ and X| are such that 

MoiXjDoas ^ t{d)t{X) => Doas = t{d)M+{X)t{X) ^ as=t{d)D^^a 
piXs) - t(Xs) - Mo{Xs)as = t(d)i(X) - t{d)MoiX)DoD^'a - t(d)p(X) 

If we consider the upper bound ([6]) of Step 2 computed for the scaled empirical 
points Xg, straightforward computations show that 

/-A/o(X5)A/+(Xs) = /-Mo(X)A/+(X) 
dkt{Xs) - M9od^s)as = ^k.t(X)-Afao.(X)« 

It follows that t satisfies condition fB]) with input set X^ if and only if 

n 4- 

|t(d)||p(X)| > \t{d)\ |/-A/o(X)A/+(X)|^-^ 9fct(X)-A/a.o(X)a 



fe=i 



141 



that is if and only if t satisfies condition ^ with input set X^ since Sk = \dk\ek- 
We conclude that the NBM algorithm puts t into O processing the input X^ if 
and only if t is added to O using the input X|. 



Let us consider the set Xt of the translated points. 



Given a term r — 



c^", there exist (see [15]) a set i? = {r 



\r} of 



power products and a set {7^ : 7^ = 7j(ui, . . . , w„)} of coefficients such that for 
each p = (ci, . . . , c„) and p = (ci +vi,.. . ,c„ + v„) 

r{p) = r(p)+ ^ IjTjip) 



Furthermore, let Fi 



(vi 



,v„) ■■ 



rjGR 

be a function such that 



TjeR 



J -^nj 



Since i^(i,i,...,i,„)(p) = for each point p G M" we obtain 







dF, 



(Vi,...,Vn) 



dxk 
that is, using our notation. 



ip) 



dr 
dxk 



dr 
dxk 



®-£-w-E7.£^(p) 



r,e-R 



' dxk ' 



dkr{p) =dkr{p) + ^ •jjdkrjip) 



11 



Now, let us consider at the current step the set O and the power product t. By 
construction O is factor closed, so that for each r € O U C[0] the set i? is a 
subset of O. Since t e C[0], we have 

t{pi) = t{pi) + ^ >^]tj{pi) and dkt{pi) = dkt{pt) + ^ X-jdktj{pi) 

SO that, denoting by A the vector which consists of the values Xj, 

t{XT) = t(X) + Afo(X)A and dktiXr) = ^^^(X) + A/a,o(X)A 

Analogously, analyzing each column of the matrices Mo{Xt) and Mqi^O{^t)^ 
there exists a square matrix A such that 

Moi^r) - A^o(X) + Afo(X)A and A%o(Xt) = Ma,o(X) + A%o(X)A (7) 

The least squares problems Mci(X)a = i(X) and A/c)(XT)aT = ^(Xt) solved 
with the input sets X^ and XJ. are such that 

Mo(X)(/ + A)aT=i(X)+Mo(X)A ^ (/ + A)aT = « + A 

^t) = i(XT) - A/o(XT)aT = t(X) + A/o(X)A - Afo(X)(a + A) = p{X) 



Since the residual of least squares problem is invariant w.r.t. the translation and 
Moi^) has full rank then Mo{Xt) is a full rank matrix too. It follows from ([7]) 
that Mo pi) {I + A) = MoO^t), and so / + A is a non singular matrix. 

If we consider the upper bound ([6]) of Step 2, computed for the translated 
empirical points X^^ straightforward calculations lead to 

/ - Mo(Kt)M^(Kt) = I- MoOC)M+(K) 

Furthermore since 

dktiXr) - MooA^tW = 5fet(X) + A%o(X)A - A/ao.(X)(/ + K)aT = 
dkt{X) + A/a,o(X)A - A/ao.(X)(a + A) = dkt{X) - Afao.(X)a 

it follows that t satisfies condition ^ with input X^ if and only if t satisfies 
condition ^ with input X"^. We conclude that the NBM algorithm puts t into 
O processing X^ if and only if t is added to O using X^. (} 

In order to analyze the stability properties of the order ideal O we recall 
some basic concepts (see [5]). 

Definition 4.2 An order ideal O is stable w.r.t. X^ if the evaluation ma- 
trix Mq (X) ha.s full rank for each admissible perturbation X o/ X"^ . 

Heuristically speaking an order ideal O can be considered stable w.r.t. the 
data uncertainty if the linear independency of the evaluation vectors of its ele- 
ments is not affected by slight perturbations of X. 
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It is well known that each order ideal is stable providing the values of e are 
sufBciently small. Nevertheless in our problem the tolerance e is given a priori 
and then not all the order ideals turn out to be stable. 

By the very nature of the NBM algorithm, no formal results about the 
stability of O can be stated. In fact, when the numerical independence of 
{r(X) : r E O} is tested using condition ([5]), Theorem 13.51 ensures that O is 
stable. Unfortunately, for implementative reasons, the NBM algorithm checks 
the numerical independence of {r(X) : r G O} using the first order approxima- 
tion ([6]) of ([5]). So the stability of O is not guaranteed. 

Nevertheless, in most cases, the numerical tests show that the upper bound ([6|) 
is satisfied with a wide margin, widely greater than O(e^j): then, although we 
cannot have the complete certainty, there is an high probability that the order 
ideal O is stable. 

We recall that (see [2]) it is possible to compute a stable order ideal by 
using the SOI algorithm. Its elevated computational cost, widely greater than 
the computational cost of the NMB algorithm, makes the SOI algorithm not 
particularly suitable for all who are not mainly interested in stability. 

However, the possibility of comparing the results of the SOI and the NBM 
algorithms points out a comforting behaviour of the NBM algorithm. In fact 
in several numerical tests the order ideals computed by the algorithms coin- 
cide. This result supports the fact that the NBM algorithm, although without 
certainty, often gives stable order ideals. 

5 Properties of the polynomial set Q 

First of all, we formalize the idea of almost vanishing polynomials introducing 
the following definition. Theorem 14.11 allows to restrict our attention to set of 
points whose coordinate belong to [—1, 1]. 

Definition 5.1 Given a set X*^ of empirical points whose coordinates belong to 
[—1, 1] a polynomial g, with coefficient vector c, is almost vanishing at X ij 

UM\2 



c 2 



< 0{eM) 



Obviously in the general case Q is not a basis of I(X) , since Q can contain 
polynomials that do not exactly vanish at X. However the following theorem 
shows that Q exhibits interesting properties w.r.t. the data uncertainty. 

Theorem 5.2 Let X"^ he a set of s empirical points and let X be an admissible 
perturbation o/X. The polynomial set Q satisfies the following properties. 

PI If g is a polynomial of Q of degree Aeg{g) and coefficient vector c, then 
"^(^^"'<.deg(5)5:e. and 



l^ll^ 



Therefore, g is almost vanishing at 
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ll.9(X)||2 

l|c||2 


'- <2s 


deg(g)> ]ek + 0{eli 
fc=i 


at X an 


id at X. 





P2 If the zero set of G is an admissible perturbation X such that MoilC) has 
full rank, then Q is the a-Grobner basis ofX{K). 

P3 // ^O ~ s, each polynomial g € Q corresponds to a unique polynomial 
gi, of the O -border basis o/X(X) such that the support of g is a subset of 
the support of gi, . Furthermore, if c and Cf, are respectively the coefficient 
vectors of g and g^ then 



l|cfc-[c,0...0]| 



II " 

^ < deg(g)||Afo(X)||2||Afo'WI|2 5]£fc 



Proof. 



PI Let us consider Step 3 of the NBM algorithm where the polynomial g is 
computed. Let t be the monomial analyzed at such step and let Ot be the 
order ideal obtained at the previous ones. Since the polynomial g is added 
to G if the residual p(X) of the least squares problem A/cij(X)a — i(X) 
does not satisfy condition ((6]) and since g^K) = p(X) we have 

n 

||5(X)||2 < ||/-Mo,(X)M+ (X)||2^£fc||afef(X)-A%c),(X)a||2 

fc=i 

First of aU we prove that ||/ - Mc>,(X)M(J,^(X)||2 = 1. In fact let A be 
a, p X q full rank matrix A, p > q and let A — t/SV^* be its singular 
values decomposition. It is well known (see [B]) that U and V are square 
orthonormal matrices and E is a block matrix of the form S* = [Si, 0], 
where Si is the square diagonal matrix of the singular values, and so 
|l/-AA+||2 = ||/-S(E*E)-iS||2 = l. 

Later on we denote by Mfc the matrix [9fet(X) Mg^Oj(X)], which consists 
of the vectors dkr{X) with r & Ot Li {t}. Obviously, if deg(a;fe,r) = 
the corresponding column of Mk is a null vector. Moreover, for each 
g, r e Ot U {t} such that q ^ r, deg{xk,r) ^ and deg{xk,q) ^ 0, we 
have dkr/deg{xk,r) ^ dkq/deg{xk,q) and, since Ot is factor closed and 
teC[0],dkr/dixk,r)eOt. 

It follows that each column dkr{lC) of Alk is a null vector or it corresponds 
to a unique column of Mc)^{X) multiplied by deg(a;fe,r'). Since ||Mfc||2 is 
equal to the norm of its submatrix consisting of the non zero columns and 
since deg(a;fe,r) < deg(g) we have that 

||A4||2<deg(5)||Mo,(X)||2 
Finally, since c = [1, —a]* is the coefficient vector of g, we have 

\\dkt{X) - Ma^oA^hh = pfecl^ < deg(5)llMo,(X)||2||c||2 (8) 
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so that 



|ff(X)||2<deg(ff)||A/o.(X)||2||c||2^efc 



(9) 



fc=i 



Since the coordinates of the points belong to [—1, 1], we have |jMc)t(X)||2 < 
||Afc)j(X)||f^ < s, where || • \\f is the Frobenius matrix norm. Then the 
first upper bound of PI foUows immediately. 

Further, in order to show the result about (7(X), note that 



So from Lemma 



g{X) = g{X) + At- AM a 
and Lemma 13.41 we have 



= g(X) + J2Ek (5fct(X) - Ma,o, (X)a) + 0(6 

k=l 

and, computing the norm of .g(X), we obtain 

n 

\\g{X)h < ||ff(X)||2 + deg(5)||Afo.(X)||2||c||^e. + 0(£ 



m) 



fe=i 



< 



2 s deg(5)||c||2;^efe + 0(ei,) 



k=l 



that is the second upper bound of PI. 

P2 If the zero set of G is an admissible perturbation X, since the residuals 
associated to the elements of G vanish at X and Mq (X) has full rank, the 
NBM algorithm computes the polynomial set G with input set X and tol- 
erance e = (0, . . . , 0). Then Property P2 follows immediately because the 
NBM algorithm with a zero tolerance coincides with the exact Buchberger- 
MoUer one. 

P3 Since ^O = s then O is the quotient basis of P/X(X) and so there exists 
the O-border basis of X(X) (see [TIT). By construction, each polynomial 
g € G with leading term t and support contained in {t}UOt corresponds to 
a polynomial gb of the O-border basis of I(X) whose support is contained 
in {t} U O. If we order the elements of Ot and O in an increasing way 
w.r.t. a, then the columns of Mo^{X) coincide with the first ifOt columns 
of AlaCK) and the coefficient vectors c = [1, —a] of 5 and Cg = [1, — /3] of 
gb obey \\cg — [c, . . . , 0]|| = ||/3— [a, . . . 0]||. Moreover they are such that 



a 




= f (X) + p(X) 



Mo(X)/3 = t(X) 



Afo(X) 
Then we obtain 

Q - (3 = Ma\X)piX) 
From p(X) = .g(X) and the upper bound ([5]), the thesis follows. 



♦ 
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Let us point out two pleasant properties of the polynomial set G easily fol- 
lowing by PI and P3. 

Property PI immediately implies that G can contain almost vanishing poly- 
nomials even if the coordinates of X do not belong to [—1, 1]. In fact it is 
sufficient that the coordinates of X are "not too elevated" (see the examples of 
Section 6). 

In the case when the condition number |jM(^^(X)||2||Mc)(X)||2 (see [6]) of 
the matrix Mci(X) is "not too elevated", Property P3 implies that g is "close" 
to a polynomial gb vanishing at X. Then X can be considered a pseudozero 
set oi G, in the sense given by Stetter (see [HI [13]). 

The formal results proved above allow us to justify in details the heuristic 
properties of G described in the Introduction and in particular the reasons why 
G characterizes the input points X. 

First of all note that Property PI implies that each element g oi G assumes 
small values, and then it is almost vanishing, at X and at each admissible per- 
turbation (of course w.r.t. the norm of its coefHcient vector). 

Moreover, we recall that, by construction, each element g oi G with leading 
term t corresponds to a least squares problem Mcij(X)a = i(X), Ot C O, 
whose residual p does not satisfy condition ([Hj). Note that, since Theorem 13.51 
involves only sufficient conditions on the residual, the fact that p does not satisfy 
condition ([Sj) gives essentially the same information of the fact that p does not 
satisfy condition ([5]). 

Given g oi G, let us suppose that there exists an admissible perturbation Xg 
such that /3(Xg) is a null vector. This is a possible case because condition ([5]) is 
not satisfied. Moreover, let us suppose that the order ideal O is stable, so that 
the matrix MotQ^g) has full rank. Then there exists a polynomial 'g, given by 
the solution of Mc)^{Xg)a = i(Xg), having the following properties: 

• 5 exactly vanishes at Xg ; 
^ has the same support of g; 



• 



gandghsive "similar" coefficients, if the condition number |lAfc)(X)|l|lAfQ(X) 
is "not too elevated" . In fact from relations ^ and ^ we have 

||[1, -S] - [1, -a]h ^ j:l^,e,\\d,tiX)-M9,o{X)a\\ ^ 



ll[l, -«]l|2 - ||[1, -«]||2 

\\M+{X)\\\\Mam\deg{g)Y,ek 

fc=i 

In this sense g can selects a geometrical configuration X^ of points, close to 
X, that can be considered an "approximate" representation of the input points 
independent of the data errors. 

Furthermore, as we will show in the examples of Section 6, it can happen that 
the whole set of polynomials g of G selects a unique geometrical configuration X. 
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Therefore the polynomials g constitute a Grobner basis of X(X) . We can then 
conclude that G can be viewed as an approximation of a Grobner basis of an 
ideal of points X close to X and the set X can be considered as a possible "exact" 
configuration corresponding to the absence of data uncertainty. 

We point out that, once again as shown in Section 6, it can happen that 
each g coincides with g and then G itself is a Grobner basis for X(X) . 

Let us conclude this section with a short recall of the open problems already 
presented in the Introduction. They are essentially related to the existence and 
possibly the determination of X. The numerical examples show that there are 
cases when X does not exist. This seems to be due to two possible causes. 
One is because the NBM algorithm could not recognize a possible element of O 
so that a polynomial g which never vanishes at any admissible perturbation is 
added to G- The second reason is when the points of X are close to different 
incompatible geometrical configurations. However, in our numerical examples, 
in this case the NBM algorithm explicitly detects these incompatible geometrical 
configurations. 

6 Numerical examples 

The following numerical tests are performed using a prototype version of the 
NBM algorithm. An improved version of it will be included soon in CoCoALib [1]. 
In all the examples the term ordering DegLex, y < x is used; in addition, the 
coordinates of the points and the coefficients of the polynomials are displayed 
with a finite number of digits, but all computations are performed in exact 
arithmetic using CoCoA 4.7 ^. 

In Example 16.11 the NBM algorithm computes an exact Grobner basis of a 
vanishing ideal of an admissible perturbation. 

Example 6.1 Given the same data of Example 11.11 that is the set of points 
X = {(1, 1), (3, 2), (5.1, 3)}, if the tolerance is e = (0.15, 0), the NBM algorithm 
computes the quotient basis O = {1, y, y^} and the polynomial set G- 



gi = .T-2.05y+1.06 
52 = y^ — 6y^ + lly — 6 



S — ■^ ^ .,3 R„,2 



The polynomial (72 vanishes at X while gi, with coefficient vector ci, is almost 
vanishing at X since ||gi(X)||/||ci|| = 0.0162. 

Since G is the cr-Grobner basis of X(X) which corresponds to the admissible 
perturbation X = {(0.983, 1), (3.03, 2), (5.083, 3)} consisting of ahgned points, 
we conclude that the points X are misaligned because of data inaccuracy. (} 



In Example 16. 21 a set X of 20 points close to a circumference is processed and 
the NBM algorithm detects this geometrical configuration. 
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Example 6.2 Let X be a set of points obtained varying the coordinates of a 
set of 20 points lying on the circumference C of equation x^ + y^ — 1 = 0, with 
componentwise perturbations less than 10""*. 

The cr-Grobner basis of I(X) does not detect that the points X are close to a 
circumference. On the contrary, the NBM algorithm, processing the set X with 
tolerance e = (0.0001, 0.0001), computes the stable quotient basis O 

O = {1, y, X, y^, xy, y^,xy'^, y^, xy^ , y'' , xy^ , y^ , xy^,y'^, xy^,y^,xy'^, y^, xy^,y^°} 

and the set Q of polynomials 

= a;2 + 0.99999^2 -1.00002 

= xy^ - 2.00006x2/^ + 1.31256x2;^ ~ 0.31251xy^ + 0.01953x2; 

= y" - 3.000062/9 + 3.31262/^ - 1.62512/^ + 0.33202;^ - 0.01952/ 

The set Q is not a Grobner basis, but (/2 and 173 vanish at X and gi, with 
coefficient vector ci, is almost vanishing at X since |j(7i(X)||/||ci|| w 10^''. 

Moreover, since the coefficient vector of gi are close to those of the circum- 
ference C we conclude that the elements of X arc "almost lying" on C. (} 

In Example 16.31 the NBM algorithm processes the same set of points with 
two different tolerances. In the first case it detects two incompatible geometrical 
configurations close to X. In the second case, choosing a smaller tolerance, the 
NBM algorithm computes a set Q very similar to a Grobner basis of X(Xi), 
where Xi is an admissible perturbation of X. 

Example 6.3 Given the set X of points 

X = {(1, 6), (2, 3), (2.449, 2.449), (3, 2), (6, 1)} 

we consider two different tolerances. 

Firstly if e = (0.018, 0.018) the NBM algorithm computes the stable quotient 
basis O — {l,y,x,y^,y^} and the set G of polynomial which is not a basis of a 
vanishing ideal since its zero set is empty: 

= xy + 0.000082/2 - 0.00064a; - 0.001252/ - 5.99501 

= a;2 + 0.991992/2 - 11.94095a: - 11.885502/ + 46.54436 

= 2/"* - 14.4772/^ + 76.72412/2 - 14.8620a; - 188.41942/ + 214.3446 

In this case gi and 172 highlight that the points of X almost lye on the hyperbola 
xy — 6 and on the circumference x"^ +y'^ — \2x — 12y + 47. In fact we have that 
both sets of points 

Xi = {(1,6),(2,3),(V6,%/6),(3,2),(6,1)} and 

X2 - {(1,6), (2, 3), (6 -2.5^2,6-2.572), (3, 2), (6,1)} 

are admissible perturbations of X. Nevertheless the configurations correspond- 
ing to Xi and X2 are incompatible, since #X = 5 while the intersection between 
an hyperbola and a circumference consists of at most 4 points. 
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If we choose a smaller tolerance, the configuration of points near to the 
circumference is not detected by the algorithm. In fact, if £ = (0.001, 0.001) we 
obtain the stable quotient basis O = {1, y, x, y^, x^} and the set G of polynomial, 
with an empty zero set: 



51 = 


= xy + O.OOOOSy^ 


~ 0.00064X - 0.001252; - 5.9950 




52 = 


= y3 _ 2.3444a;2 - 


14.3444y2 + 34.1336a; + 75.13367; - 


-182.1901 


53 = 


= x^ - 14.3444x2 


- 2. 3444^2 + 75.1336a; + 34.1336?/ - 


- 182.1901 



In this case the c-Grobner basis GBi of X(Xi) 

{xy~6 
j/3 - 2.4494a;2 - U.UMy^ + 35.3938a; + 76.3938y - 187.1260 
x^ - 14.4494a;2 - 2.4494y2 + 76.3938a; + 35.3938y - 187.1260 

consists of polynomials "similar" to the elements of Q. Since Xi is an admissible 
perturbation also w.r.t. the tolerance (0.01, 0.01) then Q highlight that the points 
X almost lye on a hyperbola. <D 



Example 16.41 shows that the term ordering a is only a computational tool 
for building a factor closed set. In fact, given the set X, the NBM algorithm 
computes the order ideal O which cannot be obtained by the exact Buchberger- 
MoUer algorithm working on X with any term ordering. 

Example 6.4 Let X = {(1.1,1.1), (0.9,-1.1), (-0.9,0.9), (-1.1,-0.9)} be 
the input points and let e — (0.12, 0.12) be the tolerance. Since the vector space 
P/X(X) has dimension 4, the possible quotient bases are 

Oi = {l,a;,a;2,a;3} O2 ^ {l^y^^y''} 
©3 = {1, 2/, X, a;2} O4 = {1, y, x, y^} O5 = {1, y, a;, xy} 

Each quotient basis Oj, j — I .. .4, is associated to the t7j-Gr6bner basis of X(X), 
where ai = a2 = Lex with y > x or x > y respectively, and a^ — ai — DegLex 
with y > X 01 X > y respectively. Nevertheless these sets are not stable quotient 
bases, since each evaluation matrix Ma-{^), j — 1...4, is singular for the 
admissible perturbation X = {(1, 1), (1,-1), (—1,1), (-1,-1)}. 

Vice versa O5, computed by the NBM algorithm, cannot be obtained using 
the exact Buchberger-MoUer algorithm w.r.t. any term ordering r. In fact the 
vector t(X), with t = x"^ or t = y"^ is independent of {r(X) : r e {l,a;,y}} so 
that O3, if a;2 <r xy, or O4, if y^ <^ xy, is built. It follows that the set G 
computed by the NBM algorithm 

y^ - 0.19998a; + 0.01980y - 1.01 

x^ - 0.20199a;?; + 0.00201a; + 0.01999y - 0.98980, 



S ~ -) „2 



is not the r-Grobner basis of X(X), for any term ordering r. Nevertheless, G 
is the cr-Grobner basis of I(X), where the zero set X of ^ is the admissible 
perturbation: 

X = {(1.099, 1.099), (0.899, -1.100), (-0.899, 0.901), (-1.099, -0.898)} <> 
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